Mini pruebas:
1 ciclo de asimilación dentro del experimento 20181122, para el 20181123000000 donde hay disponible observaciones del AMSU-A en NOAA15 y NOAA 18.
- clearsky: configuración de GSI por defecto, solo asimila radianzas en cielo claro
- allsky: activo la opción allsky modificando las variables de entrada para el CRTM en el archivo anainfo (CRTM use = 12 en vez de 10).
- noQC6: desactivo control de calidad asociado al scattering index (resta entre CH1 y CH15 para agua, valor fijo para tierra según configuración inicial).
- QC6ecmwf: modifico el cálculo del SI sobre tierra siguiendo la nota técnica del ECMWF (Weston et.al. 2019).
- QC6ecmwfall: como QC6ecmwf pero para allsky
- QC4ori: configuración por defecto para factch4 (< 0.5 no rechaza), esto es un umbral para la diferencia entre la observación y el guess del CH4 con algunas constantes en el medio.
- QC4new: se relaja el umbral a factch4 < 1.
files <- list.files("analisis/radianzas", pattern = "amsua_n15", full.names = TRUE) %>%
.[!str_detect(., "conv")]
diag <- map(files, function(f){
meta <- unglue::unglue(f, "analisis/radianzas/asim_{sensor}_{plat}_{date}.ensmean_{exp}")
# print(f)
out <- fread(f)
# .[V10 == 1] %>%
if (file.size(f) != 0) {
out[, date := ymd_hms(meta[[1]][["date"]])] %>%
.[, exp := meta[[1]][["exp"]]]
}
out
}) %>%
rbindlist()
colnames(diag) <- c("sensor", "channel", "lat", "lon", "press", "dhr", "tb_obs", "tbc", "tbcnob",
"errinv", "qc", "emis", "tlapchn", "rzen", "razi", "rlnd", "rice", "rsnw", "rcld",
"rcldp", paste0("pred", seq(8)), "date", "exp")
diag[qc == 0, .N, by = .(channel, exp)] %>%
ggplot(aes(factor(channel), exp)) +
geom_raster(aes(fill = N)) +
scale_fill_viridis_c() +
labs(x = "channel", y = "", subtitle = "QC = 0") +
theme_minimal()

diag[qc == 0 & errinv !=0, .N, by = .(channel, exp)] %>%
dcast(exp ~ channel) %>%
knitr::kable(caption = "Cantidad de observaciones potencialmente asimilables (QC == = y errinv != 0) para cada prueba y cada canal del AMSU-A")
## Using 'N' as value column. Use 'value.var' to override
Cantidad de observaciones potencialmente asimilables (QC == = y errinv != 0) para cada prueba y cada canal del AMSU-A
| QC4new |
127 |
129 |
129 |
129 |
129 |
404 |
759 |
420 |
126 |
| QC4old |
105 |
107 |
107 |
107 |
107 |
404 |
759 |
421 |
104 |
| QC4ori |
103 |
105 |
105 |
105 |
105 |
404 |
759 |
420 |
102 |
| QC6ecmwf |
99 |
101 |
101 |
101 |
101 |
403 |
709 |
410 |
98 |
| QC6ecmwfall |
184 |
183 |
176 |
186 |
192 |
402 |
716 |
413 |
181 |
| allsky |
147 |
146 |
139 |
149 |
155 |
328 |
689 |
396 |
144 |
| clearsky |
65 |
67 |
67 |
67 |
67 |
328 |
679 |
389 |
64 |
| noQC6 |
99 |
101 |
101 |
101 |
101 |
496 |
709 |
410 |
98 |
diag[channel %in% c(1:8)] %>%
dcast(lon + lat + channel ~ exp, value.var = "qc") %>%
.[, diff := allsky - clearsky] %>%
.[diff != 0] %>%
ggplot(aes(ConvertLongitude(lon), lat)) +
geom_point(color = "darkorange", size = 0.7) +
# geom_point(aes(color = errinv)) +
scale_color_viridis_c() +
geom_mapa() +
coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
facet_grid(.~ channel) +
labs(x = "lon",
subtitle = "Ubicación de observaciones extras en allsky respecto de clearsky") +
theme_minimal()

diag[channel %in% c(4)] %>%
ggplot(aes(ConvertLongitude(lon), lat)) +
geom_point(aes(color = factor(abs(qc))), size = 0.7) +
scale_color_viridis_d() +
geom_mapa() +
coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
# facet_grid(exp~channel) +
facet_wrap(~exp, ncol = 4) +
labs(subtitle = "Control de calidad para el canal 4",
x = "", y = "",
color = "QC") +
theme_minimal()

diag[channel %in% c(4:8) & qc == 0 & errinv != 0] %>%
ggplot(aes(ConvertLongitude(lon), lat)) +
geom_point(aes(color = tbc), size = 0.7) +
scale_color_viridis_c(option = "A") +
# scale_color_paletteer_c("ggthemes::Red-Blue Diverging") +
geom_mapa() +
coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
facet_grid(channel~exp) +
labs(x = "", y = "",
color = "O-B") +
theme_minimal()

# diag[channel %in% c(1:8)] %>%
# ggplot(aes(ConvertLongitude(lon), lat)) +
# geom_point(aes(color = factor(qc))) +
# scale_color_viridis_d() +
# geom_mapa() +
# coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
# facet_grid(exp~channel) +
# theme_minimal()
# diag[channel %in% c(1:15) & errinv != 0, .N, by = .(channel, exp)] %>%
# dcast(exp~channel)
diag[, qc := abs(qc)] %>%
.[channel %in% (5), .N, by = .(exp, qc)] %>%
dcast(exp~qc) %>%
knitr::kable(caption = "CategorÃas de control de calidad por experimento para el canal 5.")
## Using 'N' as value column. Use 'value.var' to override
CategorÃas de control de calidad por experimento para el canal 5.
| QC4new |
129 |
NA |
25 |
228 |
308 |
92 |
| QC4old |
107 |
NA |
23 |
138 |
307 |
207 |
| QC4ori |
105 |
NA |
21 |
139 |
308 |
209 |
| QC6ecmwf |
101 |
NA |
18 |
142 |
310 |
211 |
| QC6ecmwfall |
192 |
3 |
21 |
87 |
311 |
168 |
| allsky |
155 |
3 |
4 |
38 |
450 |
132 |
| clearsky |
67 |
NA |
3 |
91 |
450 |
171 |
| noQC6 |
101 |
NA |
34 |
190 |
NA |
457 |
diag[channel %in% c(1:15) & qc == 0 & errinv > 0.000031623,
.N, by = .(exp)] %>%
dcast(exp~.) %>%
.[order(-.)] %>%
knitr::kable(caption = "Cantidad total de observaciones potencialmente asimilables")
## Using 'N' as value column. Use 'value.var' to override
Cantidad total de observaciones potencialmente asimilables
| QC6ecmwfall |
2633 |
| QC4new |
2352 |
| allsky |
2293 |
| QC4old |
2221 |
| noQC6 |
2216 |
| QC4ori |
2208 |
| QC6ecmwf |
2123 |
| clearsky |
1793 |
diag[channel %in% c(4:8)] %>%
ggplot(aes(ConvertLongitude(lon), lat)) +
geom_point(aes(color = errinv)) +
scale_color_viridis_c() +
geom_mapa() +
coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
facet_grid(channel~exp) +
labs(x = "", y = "") +
theme_minimal()

files <- list.files("analisis/radianzas", pattern = "analysis", full.names = TRUE) %>%
.[!str_detect(., "asim")]
update <- map(files, function(f) {
descriptores <- unglue::unglue(f, c("analisis/radianzas/analysis.ensmean_{exp}"))
ana <- ReadNetCDF(f, vars = c(p = "P", "PB", t = "T", qv = "QVAPOR",
lon = "XLONG", lat = "XLAT")) %>%
.[, p := p + PB] %>%
.[, t := tk(t, p)] %>%
.[, rh := rh(qv, p, t)] %>%
.[, td := td(qv, p) + 273.15] %>%
.[, ":="(Time = NULL,
# west_east = NULL,
# south_north = NULL,
qv = NULL,
PB = NULL)] %>%
.[, c("u", "v") := uvmet(f)] %>%
.[, ":="(exp = descriptores[[1]][["exp"]])] %>%
.[]
}) %>%
rbindlist()
guess <- ReadNetCDF("analisis/radianzas/wrfarw.ensmean",
vars = c(p = "P", "PB", t = "T", qv = "QVAPOR",
lon = "XLONG", lat = "XLAT")) %>%
.[, p := p + PB] %>%
.[, t := tk(t, p)] %>%
.[, rh := rh(qv, p, t)] %>%
.[, td := td(qv, p) + 273.15] %>%
.[, ":="(Time = NULL,
# west_east = NULL,
# south_north = NULL,
qv = NULL,
PB = NULL)] %>%
.[, c("u", "v") := uvmet("analisis/radianzas/wrfarw.ensmean")] %>%
# .[, ":="(exp = "guess")] %>%
.[]
update <- guess[update, on = c("lon", "lat", "bottom_top", "south_north", "west_east")]
update[, .(mean_ana = mean(i.t, na.rm = TRUE),
mean_guess = mean(t, na.rm = TRUE)), by = .(bottom_top, west_east, exp)] %>%
ggplot(aes(west_east, bottom_top)) +
geom_contour_fill(aes(z = (mean_ana - mean_guess))) +
scale_fill_divergent() +
labs(fill = "update T") +
facet_wrap(~ exp) +
theme_minimal()

update[, .(mean_ana = mean(i.t, na.rm = TRUE),
mean_guess = mean(t, na.rm = TRUE)), by = .(bottom_top, south_north, exp)] %>%
ggplot(aes(south_north, bottom_top)) +
geom_contour_fill(aes(z = (mean_ana - mean_guess))) +
scale_fill_divergent() +
labs(fill = "update T") +
facet_wrap(~ exp) +
theme_minimal()

update %>%
dcast(bottom_top + south_north + west_east ~ exp, value.var = "i.t") %>%
.[, diff := allsky - QC6ecmwf] %>%
.[diff != 0]
## Empty data.table (0 rows and 7 cols): bottom_top,south_north,west_east,QC6ecmwf,allsky,clearsky...